##########    ESS partner analysis
##########    Jonne Kamphorst
##########    
##########3


##### Load libraries and data ####
library(haven)
library(sjPlot)
library(ggplot2)
library(extrafont)
library(ggeffects)
library(gridExtra)
library(ggpubr)
library(stargazer)
library(dplyr)
library(tidyverse)
library(stringr)
library(viridis)
library(xtable)
library(interplot)
library(sjPlot)
library(texreg)
library(estimatr)
library(ggeffects)
library(viridis)



#font_import()
loadfonts(quiet = T, device = "pdf")
windowsFonts(Georgia = windowsFont("Georgia")) #load fonts for windows machines


## Set ggplot theme with the Georgia font
theme_set(theme_light(base_family = "Georgia",  base_size = 12))


# Set working directory and load data
#JONNE
ESS <- readRDS("ESS partner/ESS_complete.rds")

## Keep countries we are interested in
ESS <- subset(ESS, cntry == "BE" | cntry == "GB" | 
                cntry == "FR" | cntry == "NL" | 
                cntry == "IR"| cntry == "DE"  | cntry == "DK" | 
                cntry == "NO" | cntry == "SE" | cntry == "FI" | cntry == "CH" |
                cntry == "ES"  | cntry == "PT")



## Remove people who are too young to be out of university already
ESS <- subset(ESS, agea > 26)






############################### Clean and recode variables #################################

# Add generations
ESS$generation <- NA
ESS$generation[ESS$yrbrn < 1946] <- 1
ESS$generation[ESS$yrbrn > 1945 & ESS$yrbrn < 1956] <- 2
ESS$generation[ESS$yrbrn > 1955 & ESS$yrbrn < 1966] <- 3
ESS$generation[ESS$yrbrn > 1965 & ESS$yrbrn < 1976] <- 4
ESS$generation[ESS$yrbrn > 1975 & ESS$yrbrn < 1986] <- 5
ESS$generation[ESS$yrbrn > 1985 & ESS$yrbrn < 1996] <- 6





######## Clean Education variables
## Education respondent
# edulvla Harmonised variable (fewer categories-more rounds)
# edulvlb (more categories, only from ESS5)
ESS$education <- ESS$edulvla
ESS$education[ESS$education == 55 | ESS$education == 77 | ESS$education == 88 | ESS$education == 99 | ESS$education == 0] <- NA #these categories are things like 'didn't answer'. The 0 is that it can't be put into this harmonized scheme that they have for all rounds


## Education partner
# edulvlpb education partner
ESS$education_par <- ESS$edulvlpa
ESS$education_par[ESS$education_par == 55 | ESS$education_par == 77 | ESS$education_par == 88 | ESS$education_par == 99 | ESS$education_par == 0 ] <- NA #66 category not among NA(not applicable)


## Education father
# edulvlfa
ESS$education_fat <- ESS$edulvlfa
ESS$education_fat[ESS$education_fat == 55 | ESS$education_fat == 77 | ESS$education_fat == 88 | ESS$education_fat == 99 | ESS$education_fat == 0] <- NA


## Education mother
# edulvlma
ESS$education_mot <- ESS$edulvlma
ESS$education_mot[ESS$education_mot == 55 | ESS$education_mot == 77 | ESS$education_mot == 88 | ESS$education_mot == 99 | ESS$education_mot == 0] <- NA


########## Recode education variables
# Add dummies for education
# Education respondent
ESS$education <- as.numeric(ESS$education)
table(ESS$education)
ESS$education_university_dummy <- ifelse(ESS$education > 4, 1, 0) # 1 is university 0 is not university

# education_par partner
ESS$education_par <- as.numeric(ESS$education_par)
ESS$education_par_university_dummy <- ifelse(ESS$education_par > 4, 1, 0)

# education_fat  father
ESS$education_fat <- as.numeric(ESS$education_fat)
ESS$education_university_father_dummy <- ifelse(ESS$education_fat > 4, 1, 0)

# education_mot  mother
ESS$education_mot <- as.numeric(ESS$education_mot)
ESS$education_university_mother_dummy <- ifelse(ESS$education_mot > 4, 1, 0)
sum(table(ESS$education_mot)) == sum(table(ESS$education_university_mother_dummy)) # should be true

# create a dummy variable for both parents being at university
ESS$education_university_parents_dummy <- ifelse(ESS$education_university_mother_dummy == 1 & 
                                                   ESS$education_university_father_dummy == 1,  1, 0) 

#Remove cases where we don't have the info for parents
ESS$education_university_parents_dummy[is.na(ESS$education_university_mother_dummy)] <- NA
ESS$education_university_parents_dummy[is.na(ESS$education_university_father_dummy)] <- NA





######### Add the main instrument for educational sorting
ESS$sorting <- NA
ESS$sorting[ESS$education_university_dummy == 1 & ESS$education_university_parents_dummy == 1 & ESS$education_par_university_dummy == 1 | 
              ESS$education_university_dummy == 1 & ESS$education_university_parents_dummy == 1 & ESS$partner == 2] <- "Pure High"
# R uni, M uni, F uni, P uni ;; R uni, M uni, F uni, No P = pure uni  

ESS$sorting[ESS$education_university_dummy == 1 & ESS$education_university_parents_dummy == 1 & ESS$education_par_university_dummy == 0 | 
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 & ESS$education_par_university_dummy == 1|
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1 & ESS$education_par_university_dummy == 1|
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1 & ESS$partner == 2|
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 & ESS$partner == 2] <- "High weak ties"              
# R uni, M uni, F uni, P low ;; R uni, M uni, F low, P uni ;; R uni, M low, F uni, P uni ;; R uni, M low, F, high, no P ;; R uni, M high, F low, no P = "high weak ties"

ESS$sorting[ESS$education_university_dummy == 1 & ESS$education_university_parents_dummy == 0 & ESS$education_par_university_dummy == 0 | 
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 & ESS$education_par_university_dummy == 0|
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1 & ESS$education_par_university_dummy == 0|
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 0 & ESS$education_par_university_dummy == 1|
              ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 0 & ESS$partner == 2] <- "High strong ties"              

# R uni, MFP low ;; R uni, M uni, F low P low ;; R uni, M low, F uni P low; R uni, M low, F low P uni ;; R uni, MF low, no P = High strong ties

ESS$sorting[ESS$education_university_dummy == 0 & ESS$education_university_parents_dummy == 0 & ESS$education_par_university_dummy == 0 | 
              ESS$education_university_dummy == 0 & ESS$education_university_parents_dummy == 0 & ESS$partner == 2] <- "Pure Low"

ESS$sorting[ESS$education_university_dummy == 0 & ESS$education_university_parents_dummy == 0 & ESS$education_par_university_dummy == 1 | 
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 & ESS$education_par_university_dummy == 0|
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1 & ESS$education_par_university_dummy == 0|
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1 & ESS$partner == 2|
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 & ESS$partner == 2] <- "Low weak ties" 

ESS$sorting[ESS$education_university_dummy == 0 & ESS$education_university_parents_dummy == 1 & ESS$education_par_university_dummy == 1 | 
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 & ESS$education_par_university_dummy == 1|
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1 & ESS$education_par_university_dummy == 1|
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 1 & ESS$education_par_university_dummy == 0|
              ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 1 & ESS$partner == 2] <- "Low strong ties"  

ESS$sorting <- as.factor(ESS$sorting)







##### Add robustness intstrument for lower educted sorting ####
ESS$sorting_robst <- ESS$sorting
ESS$sorting_robst <- as.character(ESS$sorting_robst)
ESS$sorting_robst[ESS$sorting == "Pure Low" & (ESS$education == 1 | ESS$education == 2)] <- "Pure Low - Low"
ESS$sorting_robst[ESS$sorting == "Low weak ties" & (ESS$education == 1 | ESS$education == 2)] <- "Low weak ties - Low"
ESS$sorting_robst[ESS$sorting == "Low strong ties" & (ESS$education == 1 | ESS$education == 2)] <- "Low strong ties - Low"
ESS$sorting_robst[ESS$sorting == "Pure Low" & (ESS$education == 3 | ESS$education == 4)] <- "Pure Low - Med"
ESS$sorting_robst[ESS$sorting == "Low weak ties" & (ESS$education == 3 | ESS$education == 4)] <- "Low weak ties - Med"
ESS$sorting_robst[ESS$sorting == "Low strong ties" & (ESS$education == 3 | ESS$education == 4)] <- "Low strong ties - Med"






############# Add sorting for parents only ##############

## Add the sorting code for parents only
ESS$sorting_parents <- NA
ESS$sorting_parents[ESS$education_university_dummy == 1 & ESS$education_university_parents_dummy == 1] <- "Pure High" 
ESS$sorting_parents[ESS$education_university_dummy == 1 & ESS$education_university_parents_dummy == 0] <- "High strong ties" 
ESS$sorting_parents[ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 |
                      ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1] <- "High weak ties" 

ESS$sorting_parents[ESS$education_university_dummy == 0 & ESS$education_university_parents_dummy == 0] <- "Pure Low" 
ESS$sorting_parents[ESS$education_university_dummy == 0 & ESS$education_university_parents_dummy == 1] <- "Low strong ties" 
ESS$sorting_parents[ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 |
                      ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1] <- "Low weak ties" 

ESS$sorting_parents <- factor(ESS$sorting_parents, levels=c("Pure High", "High weak ties", "High strong ties",
                                                             "Low strong ties", "Low weak ties", "Pure Low"))     



## Sorting version for the plots
ESS$sorting_parents1 <- NA

ESS$sorting_parents1[ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 1 | 
                       ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 |
                       ESS$education_university_dummy == 1 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1] <- "High res High parent(s)" 

ESS$sorting_parents1[ESS$education_university_dummy == 1 & ESS$education_university_father_dummy == 0 & ESS$education_university_mother_dummy == 0] <- "High res Low parents" 


ESS$sorting_parents1[ESS$education_university_dummy == 0 & ESS$education_university_father_dummy == 0 & ESS$education_university_mother_dummy == 0] <- "Low res Low parents" 

ESS$sorting_parents1[ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 1 |
                       ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 1 & ESS$education_university_father_dummy == 0 |
                       ESS$education_university_dummy == 0 & ESS$education_university_mother_dummy == 0 & ESS$education_university_father_dummy == 1 ] <- "Low res High parent(s)" 

ESS$sorting_parents1 <- factor(ESS$sorting_parents1, levels=c("High res High parent(s)", "High res Low parents", "Low res High parent(s)", "Low res Low parents"))  







############# Add sorting for partner only ##############
ESS$sorting_partner <- NA
ESS$sorting_partner[ESS$education_university_dummy == 1  & ESS$education_par_university_dummy == 1] <- "Pure High"
# R uni, M uni, F uni, P uni ;; R uni, M uni, F uni, No P = pure uni  

ESS$sorting_partner[ESS$education_university_dummy == 1 & ESS$education_par_university_dummy == 0] <- "High ties"              
# R uni, M uni, F uni, P low ;; R uni, M uni, F low, P uni ;; R uni, M low, F uni, P uni ;; R uni, M low, F, high, no P ;; R uni, M high, F low, no P = "high weak ties"


ESS$sorting_partner[ESS$education_university_dummy == 0 & ESS$education_par_university_dummy == 0] <- "Pure Low"

ESS$sorting_partner[ESS$education_university_dummy == 0 & ESS$education_par_university_dummy == 1] <- "Low ties" 

ESS$sorting_partner <- factor(ESS$sorting_partner, levels=c("Pure High", "High ties", "Low ties", "Pure Low"))






############# Add sorting for partner w singles only ##############
ESS$sorting_partner_wsingle <- NA
ESS$sorting_partner_wsingle[ESS$education_university_dummy == 1  & ESS$education_par_university_dummy == 1] <- "Pure High w par"
# R uni, M uni, F uni, P uni ;; R uni, M uni, F uni, No P = pure uni  

ESS$sorting_partner_wsingle[ESS$education_university_dummy == 1 & ESS$education_par_university_dummy == 0] <- "High ties w par"              
# R uni, M uni, F uni, P low ;; R uni, M uni, F low, P uni ;; R uni, M low, F uni, P uni ;; R uni, M low, F, high, no P ;; R uni, M high, F low, no P = "high weak ties"


ESS$sorting_partner_wsingle[ESS$education_university_dummy == 0 & ESS$education_par_university_dummy == 0] <- "Pure Low w par"

ESS$sorting_partner_wsingle[ESS$education_university_dummy == 0 & ESS$education_par_university_dummy == 1] <- "Low ties w par" 

ESS$sorting_partner_wsingle[ESS$education_university_dummy == 1 & ESS$partner == 2] <- "Pure High no par" 

ESS$sorting_partner_wsingle[ESS$education_university_dummy == 0 & ESS$partner == 2] <- "Pure Low no par" 



ESS$sorting_partner_wsingle <- factor(ESS$sorting_partner_wsingle, levels=c("Pure High w par", 
                                                                            "Pure High no par",
                                                                            "High ties w par",
                                                                            "Low ties w par",
                                                                            "Pure Low no par",
                                                                            "Pure Low w par"))







#################   Outcome variables    ###########################
### coded as (dummies): 1 is progressive, 0 is conservative


### Migration variables 
## These variables are based on Cavaille and Marshall. Cavaille discusses three different ways to code these variables. 
## We have so far coded everything as 1 being progressive and 0 being consevative. Cavaille does this the other way around
## Of course, this does not influence the categories. 
## So that these also capture being anti-immigrant. This also fits more with what the literature tries to explain. 
##     Cavaille's coding:
##  (1) simply as dummies
##  (2) As a single dummy if you express anti-behavior on variables 1, 2, or 3 in Table 2.
##  (3) All six questions added up together and standardized
##  (4) As a fourth method she uses factor analysis. This does not make much sense in my opinion


#### Cavaille 1, 2, and 3 in Table 2. 
#  Note that cavaille uses coding with either only 'none' or also 'a few' as conservative. She said to prefer the latter in the Appendix, I follow this.
#Allow many/few immigrants from the same ethnicity as in the country
ESS$imsmetn[ESS$imsmetn == 7 | ESS$imsmetn == 8 | ESS$imsmetn == 9 ] <- NA
ESS$imsmetn_dummy <- ifelse(ESS$imsmetn > 2, 0, 1)

#Allow many/few immigrants of different race/ethnic group from majority
ESS$imdfetn[ESS$imdfetn == 7 | ESS$imdfetn == 8 | ESS$imdfetn == 9 ] <- NA
ESS$imdfetn_dummy <- ifelse(ESS$imdfetn > 2, 0, 1) # greater than 2 is the 'none or a few' category. 1 thus means that people are anti immigration

ESS$imdfetn_pro <- max(ESS$imsmetn, na.rm=T) - ESS$imsmetn



#Allow many/few immigrants from poorer countries outside Europe
ESS$impcntr[ESS$impcntr == 7 | ESS$impcntr == 8 | ESS$impcntr == 9 ] <- NA
ESS$impcntr_dummy <- ifelse(ESS$impcntr > 2, 0, 1)


### 4, 5, and 6 in cavaille table 2
# cut-off point for cavaille is always 4. 
# Immigration bad or good for country's economy
ESS$imbgeco[ESS$imbgeco == 77 | ESS$imbgeco == 88 | ESS$imbgeco == 99 ] <- NA
ESS$imbgeco_dummy <- ifelse(ESS$imbgeco > 4, 1, 0)

# Country's cultural life undermined or enriched by immigrants
ESS$imueclt[ESS$imueclt == 77 | ESS$imueclt == 88 | ESS$imueclt == 99 ] <- NA
ESS$imueclt_dummy <- ifelse(ESS$imueclt > 4, 1, 0)

# is the country made a worse or a better place to live 
ESS$imwbcnt[ESS$imwbcnt == 77 | ESS$imwbcnt == 88 | ESS$imwbcnt == 99] <- NA
ESS$immigration_dummy <- ifelse(ESS$imwbcnt > 4, 1, 0)



# everyone needs to share the same customs and traditions. Waves 1 and 7 only
ESS$pplstrd[ESS$pplstrd == 9 | ESS$pplstrd == 8 | ESS$pplstrd == 7] <- NA
ESS$pplstrd_dummy <- ifelse(ESS$pplstrd  > 3 , 1, 0) # 1"progressive"



## Immigration according to Cavaille and Marshall's Appendix. Ideal way to code variables 1-3 in paper table 2 together.
# imsmetn; imdfetn; impcntr
# if someone is *pro* on any of these three types of migrants, then you become 0 in the measure. (thus anti)
# people who are coded as 1 thus do not want one type of migrant
# Note that we can't simply change the direction of this variable as Cavaille uses it because she is worried about
# social desirability bias where people try to be 'good' by allowing one type of migrant
ESS$immi_anti_cavaille <- ifelse(ESS$imsmetn_dummy == 0 | ESS$imdfetn_dummy == 0 |
                              ESS$impcntr_dummy == 0, 1, 0)

ESS$immi_pro_cavaille <- 1- ESS$immi_anti_cavaille





        
## Immigration with all variables added together as described on page 257 of paper and in Appendix A.2.
#  her 'lumped together' measure
# it's an additive scale so all other dummy variables are added together. The measures are then standardized.

# First create the lumped together variable
ESS$immi_anti_cavaille_lump <- ESS$imsmetn_dummy + ESS$imdfetn_dummy + ESS$impcntr_dummy + # variables 1,2, and 3 from Table 2
  ESS$imbgeco_dummy + ESS$imueclt_dummy + ESS$immigration_dummy # variables 4, 5, and 6

#then standardize between and within
ESS <- ESS %>%  group_by(cntry) %>% #group by country to standardize within countries
  mutate(immi_anti_cavaille_lump_sd = sd(immi_anti_cavaille_lump, na.rm=T), #calculate mean and SD for standardization
                      immi_anti_cavaille_lump_mean = mean(immi_anti_cavaille_lump, na.rm=T)) %>%
  mutate(immi_anti_cavaille_lump_standardizedwithin = #standardize within countries
           (immi_anti_cavaille_lump - immi_anti_cavaille_lump_mean)/ immi_anti_cavaille_lump_sd) %>%
  ungroup() %>%
  mutate(immi_anti_cavaille_lump_sd = sd(immi_anti_cavaille_lump, na.rm=T), #overwrite mean and SD 
         immi_anti_cavaille_lump_mean = mean(immi_anti_cavaille_lump, na.rm=T)) %>%
  mutate(immi_anti_cavaille_lump_standardizedbetween = #standardize between countries
           (immi_anti_cavaille_lump - immi_anti_cavaille_lump_mean)/ immi_anti_cavaille_lump_sd) 



### Other attitudinal outcome variables ###
## Note that here 1 is progressive


# European unification go further or gone too far, 0 gone too far, 10 go further
# in waves 2 through 9
ESS$euftf[ESS$euftf == 77 | ESS$euftf == 88 | ESS$euftf == 99] <- NA
ESS$euftf_dummy <- ifelse(ESS$euftf > 5, 1, 0) # 1"progressive" 

## Women's rights
ESS$mnrgtjb[ESS$mnrgtjb == 7 | ESS$mnrgtjb == 8 | ESS$mnrgtjb == 9] <- NA
ESS$mnrgtjb_dummy <- ifelse(ESS$mnrgtjb > 3, 1, 0)

## LGBTQ's rights
ESS$freehms[ESS$freehms == 7 | ESS$freehms == 8 | ESS$freehms == 9] <- NA
ESS$freehms_dummy <- ifelse(ESS$freehms > 2, 0, 1)

## Redistribution (Same as Paskov 2021 and Rueda foor first then morals)
# redistribution support (Government should reduce differences in income levels) USE
# in all waves
ESS$gincdif[ESS$gincdif == 7 | ESS$gincdif == 8 | ESS$gincdif == 9 ] <- NA
ESS$gincdif_dummy <- ifelse(ESS$gincdif > 2, 0, 1)

## Fairness beliefs
# smdfslv (	For fair society, differences in standard of living should be small, 1 agree strongly 5 disagree strongly)
#In waves 4 and 8
ESS$smdfslv[ESS$smdfslv == 9 | ESS$smdfslv == 8 | ESS$smdfslv == 7] <- NA
ESS$smdfslv_dummy <- ifelse(ESS$smdfslv == 1 | ESS$smdfslv == 2, 1, 0) #1 = "Equal societies need to be fair

## Fairness believes
#sblazy: Social benefits/services make people lazy
#In waves 4 and 8
ESS$sblazy[ESS$sblazy == 9 | ESS$sblazy == 8 | ESS$sblazy == 7] <- NA
ESS$sblazy_dummy <- ifelse(ESS$sblazy < 3, 0, 1)

## Political chances
# psppsgva (How much would you say the political system allows people like you to have a say in what the government does?) Higher more agree
# in waves 8 and 9 only
ESS$psppsgva[ESS$psppsgva == 9 | ESS$psppsgva == 8 | ESS$psppsgva == 7] <- NA
ESS$psppsgva_dummy <- ifelse(ESS$psppsgva == 5 | ESS$psppsgva == 4 | ESS$psppsgva == 3, 1, 0) #1 is have a say. 

## Left right self placement
ESS$lrscale[ESS$lrscale == 77 | ESS$lrscale == 88 | ESS$lrscale == 99] <- NA






############    Control variables   ###########
# ethnicity
ESS$ethnic_min <- NA
ESS$ethnic_min[ESS$blgetmg == 2] <- 0
ESS$ethnic_min[ESS$blgetmg == 1] <- 1

# Rural and urban
ESS$domicil[ESS$domicil > 5] <- NA

# Income
# only in waves 4 till 9
ESS$hinctnta[ESS$hinctnta == 77 | ESS$hinctnta == 88 | ESS$hinctnta == 99 ] <- NA  #hinctnta from ESS4 onward. 
ESS$hinctnta <- as.numeric(ESS$hinctnta)

# Gender
ESS$gndr <- as.factor(ESS$gndr)
#all rounds

# Age #all rounds
table(ESS$agea)

ESS$agea_sqrt <- ESS$agea ^2

# Employment status #all rounds
#table(ESS$pdwrk)


# Cntry
ESS$country_factor <- as.factor(ESS$cntry)

#ISCO codes
# Change into two digit isco code
ESS$isco_three <- substr(ESS$isco08,1,2)
ESS$isco_three_p <- substr(ESS$isco08p,1,2)

# Oesch class schema (added in the other OESCH files) These are the first files in the folder

# Factorize the class variables 
ESS$class8_factor <- as.factor(ESS$class8)
ESS$class5_factor <- as.factor(ESS$class5)

# Code-up parents class
ESS$occf14_factor <- as.factor(ESS$occf14b) #parents education
ESS$occm14_factor <- as.factor(ESS$occm14b)

ESS$father_class <- NA
ESS$father_class[ESS$occf14b == 1 | ESS$occf14b == 2] <- "Upper class"
ESS$father_class[ESS$occf14b == 3 | ESS$occf14b == 4 | ESS$occf14b == 5] <- "Middle class"
ESS$father_class[ESS$occf14b == 6 | ESS$occf14b == 7] <- "Skilled workers"
ESS$father_class[ESS$occf14b == 8 | ESS$occf14b == 9] <- "Unskilled workers"

ESS$mother_class <- NA
ESS$mother_class[ESS$occm14b == 1 | ESS$occm14b == 2] <- "Upper class"
ESS$mother_class[ESS$occm14b == 3 | ESS$occm14b == 4 | ESS$occm14b == 5] <- "Middle class"
ESS$mother_class[ESS$occm14b == 6 | ESS$occm14b == 7] <- "Skilled workers"
ESS$mother_class[ESS$occm14b == 8 | ESS$occm14b == 9] <- "Unskilled workers"



# factorize ess round variable
ESS$year <- NA
ESS$year[ESS$essround==1] <- 2002
ESS$year[ESS$essround==2] <- 2004
ESS$year[ESS$essround==3] <- 2006
ESS$year[ESS$essround==4] <- 2008
ESS$year[ESS$essround==5] <- 2010
ESS$year[ESS$essround==6] <- 2012
ESS$year[ESS$essround==7] <- 2014
ESS$year[ESS$essround==8] <- 2016
ESS$year[ESS$essround==9] <- 2018
ESS$year <- as.factor(ESS$year)

#Gender
ESS$female <- NA
ESS$female[ESS$gndr==1] <- 0
ESS$female[ESS$gndr==2] <- 1


# Political Interest
ESS$polintr <-  4 - ESS$polintr # Higher values more interested



##### Standardize the control variables  where that is needed
# Function to standardize a single variable (mean 0 SD of one)
stand <- function(x){
  var <- (x - mean(x,na.rm=T))/sd(x,na.rm=T)
  return(var)
}

ESS$hinctnta_01 <- stand(ESS$hinctnta)
ESS$agea10 <- ESS$agea / 10
ESS$domicil_01 <- stand(ESS$domicil)



### Add the migration index

### Migration outcome according to Kriesi and Hausermann in Beramendi 2015.
## Simple Additive index of six different immigration items: 
# imsmetn   imdfetn   impcntr   imbgeco   imueclt   imwbcnt

## Reverse the polarity of the items where higher values are not more progressive
ESS$imsmetn_rev <- max(ESS$imsmetn,na.rm=T) - ESS$imsmetn
ESS$imdfetn_rev <- max(ESS$imdfetn,na.rm=T) - ESS$imdfetn
ESS$impcntr_rev <- max(ESS$impcntr,na.rm=T) - ESS$impcntr

# Add together using ICW
# Function to standardize columns of a matrix
# where you designate a standardization group
# (e.g., the control group in an experiment)
# with "sgroup", a logical vector.

matStand <- function(x, sgroup = rep(TRUE, nrow(x))){
  for(j in 1:ncol(x)){
    x[,j] <- (x[,j] - mean(x[sgroup,j]))/sd(x[sgroup,j])
  }
  return(x)
}

# Function that takes in data in matrix format and returns
# (i) IC weights and (ii) ICW index.
# Weights can be incorporated using the "wgts" argument.
# The "revcols" argument takes a vector indicating which columns,
# if any, should have values reversed (that is, standardized 
# values are multiplied by -1) prior to construction of the index. 

icwIndex <- function(	xmat,
                      wgts=rep(1, nrow(xmat)),
                      revcols = NULL,
                      sgroup = rep(TRUE, nrow(xmat))){
  X <- matStand(xmat, sgroup)
  if(length(revcols)>0){
    X[,revcols] <-  -1*X[,revcols]
  }
  i.vec <- as.matrix(rep(1,ncol(xmat)))
  Sx <- cov.wt(X, wt=wgts)[[1]]
  weights <- solve(t(i.vec)%*%solve(Sx)%*%i.vec)%*%t(i.vec)%*%solve(Sx)
  index <- t(solve(t(i.vec)%*%solve(Sx)%*%i.vec)%*%t(i.vec)%*%solve(Sx)%*%t(X))
  return(list(weights = weights, index = index))
}



#### Select variables and apply functions
## Create a special version of the ESS DF with no missigness on key variables
ESS_immi <- ESS %>% drop_na(imsmetn_rev, imdfetn_rev, impcntr_rev, imbgeco, imueclt, imwbcnt)

matrix <- ESS_immi %>% dplyr::select(imsmetn_rev, imdfetn_rev, impcntr_rev, imbgeco, imueclt, imwbcnt)
matrix <- as.matrix(matrix)
matrix_icw <- icwIndex(matrix)
ESS_immi$immigration_icw_outcome <- matrix_icw$index
rm(matrix, matrix_icw)
min(ESS_immi$immigration_icw_outcome, na.rm=T)
max(ESS_immi$immigration_icw_outcome, na.rm=T)
sd(ESS_immi$immigration_icw_outcome, na.rm=T)
mean(ESS_immi$immigration_icw_outcome, na.rm=T)

# Change so it runs from 0 to 1
ESS_immi$immigration_icw_outcome <- (ESS_immi$immigration_icw_outcome - min(ESS_immi$immigration_icw_outcome, na.rm = TRUE)) / (max(ESS_immi$immigration_icw_outcome, na.rm = TRUE) - min(ESS_immi$immigration_icw_outcome, na.rm = TRUE))








######################## Analysis #######################
##### Main analysis

## Main outcome variables:
# euftf_dummy   - EU dummy  USE
# ESS_immi$immigration_icw_outcome
# vote.prog       - vote progressive    USE
# vote.rrp    - vote radical right    USE


#### For with parents
# create the models 

model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "sorting_parents + 
                      hinctnta_01 +  
                      female + 
                      agea10 + 
                      pdwrk + 
                      education +
                      domicil_01 +
                      ethnic_min +
                      year + 
                      class5_factor +
                      mother_class + 
                      father_class + 
                      country_factor"))
  model_full <- lm(formula=formula, data=data)
  return(model_full)
}



eu_full <- model_full("euftf", ESS)

immi1_full <- model_full("immigration_icw_outcome", ESS_immi) #ICW index

voteprog_full <- model_full("vote.prog", ESS)

voterrp_full <- model_full("vote.rrp", ESS)

redis_full <- model_full("gincdif", ESS)



### Table A4. 
# Model with variables for in the appendix
model <- texreg(l=list(eu_full, immi1_full, voterrp_full, voteprog_full, redis_full),
                stars=c(0.05, 0.01, 0.001),use.packages=F,
                digits=3, include.ci = FALSE, booktabs=T, fontsize="small",
                custom.model.names = c("EU", "Immigration", "Voting RRP", "Voting Green", "Redistribution"),
                omit.coef=c("(year)|(class5)|(mother)|(father)|(country)"), 
                custom.coef.names=c("Intercept", "High education One low parent", "High education Low parents", "Low education High parents", 
                                    "Low education One high parent", "Low education Low parents", "Income (standardized)", "Female", "Age (in 10 yrs)", "Employed (0-1)",
                                    "Education (continuous)", "Rural (standardized)", "Ethnic minority"),
                reorder.coef=c(2,3,4,5,6,7,8,9,10,11,12,13,1),
                caption.above=T,
                caption="ESS: Effect of parental education on attitudes and voting with controls",
                label="tab:main_ess_parents_withcontrols",
                threeparttable = TRUE,
                float.pos="htpb!",
                custom.gof.rows=list("Controls"= c("Yes", "Yes", "Yes","Yes", "Yes"),
                                     "Year FE"= c("Yes", "Yes", "Yes","Yes", "Yes"),
                                     "Country FE"= c("Yes", "Yes", "Yes","Yes", "Yes"),
                                     "Class FE"= c("Yes", "Yes", "Yes","Yes", "Yes"),
                                     "Parent's class FE FE"= c("Yes", "Yes", "Yes","Yes", "Yes")),
                custom.note="\\item %stars. The table shows that the education of someone's parents influences their attitudes and voting behavior. The reference category for the sorting variable are Higher educated respondents who's parents are both higher educated. The control variables (except education and the numer of ties) are standardized. EU is on a 10-step Likert scale and immigration is an ICW scale composed of 6 individual items (0-1). Redistribution attitudes are on a 5-step Likert scale.")

print(model, file = "ESS partner/tables/taba4_main_ess_parents_withcontrols.tex")
rm(model)






#### For with your partner
model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "sorting_partner_wsingle +
                     hinctnta_01 +
                     female +
                     agea10 + 
                     pdwrk + 
                     education +
                     domicil_01 +
                     ethnic_min +
                     polintr + 
                     class5_factor + 
                     mother_class + 
                     father_class + 
                     education_university_mother_dummy +
                     education_university_father_dummy + 
                     country_factor +
                     year"))
  model_full<- lm(formula=formula, data=data)
  return(model_full)
}


eu_full <- model_full("euftf", ESS)

immi1_full <- model_full("immigration_icw_outcome", ESS_immi) #ICW index

voteprog_full <- model_full("vote.prog", ESS)

voterrp_full <- model_full("vote.rrp", ESS)

redis_full <- model_full("gincdif", ESS)




#### Table A3
## Model with controls for in the appendix
model <- texreg(l=list(eu_full, immi1_full, voterrp_full, voteprog_full, redis_full),
                stars=c(0.05, 0.01, 0.001),use.packages=F,
                digits=3, include.ci = FALSE, booktabs=T, fontsize="small",
                custom.model.names = c("EU", "Immigration", "Voting RRP", "Voting Green", "Redistribution"),
                omit.coef=c("(year)|(class5)|(mother_class)|(father_class)|(country)"), 
                custom.coef.names=c("Intercept", "High res no par", "High res Low par",
                                    "Low res High par", "Low res no par", 
                                    "Low res Low par", "Income (standardized)", "Female", "Age (in 10 yrs)", "Employed (0-1)",
                                    "Education (continuous)", "Rural (standardized)", "Ethnic minority", "Political interest (0-3)", "Mother university", "Father university"),
                reorder.coef=c(2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1),
                label="tab:main_ess_partner_wcontrols",
                caption.above=T,
                caption="ESS: Effect of partner education on attitudes and voting with controls",
                threeparttable = TRUE,
                float.pos="htpb!",
                custom.gof.rows=list("Controls"= c("Yes", "Yes", "Yes","Yes", "Yes"),
                                     "Year FE"= c("Yes", "Yes", "Yes","Yes", "Yes"),
                                     "Country FE"= c("Yes", "Yes", "Yes","Yes", "Yes"),
                                     "Class FE"= c("Yes", "Yes", "Yes","Yes", "Yes")),
                custom.note="\\item %stars. The table shows that the education level of someone's partner influences their attitudes and voting behavior. The reference category for the sorting variable are Higher educated respondents who's parents are both higher educated. EU is on a 10-step Likert scale and immigration is an ICW scale composed of 6 individual items (0-1). Redistribution attitudes are on a 5-step Likert scale.")

print(model, file = "ESS partner/tables/taba3_main_ess_partner_wcontrols.tex")













############# For partners, for each outcome, a graph with predicted probabilities
# and for the means per education group
model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "education_university_dummy +
                     hinctnta_01 +
                     female +
                     agea10 + 
                     pdwrk + 
                     education +
                     domicil_01 +
                     polintr + 
                     ethnic_min +
                     mother_class + 
                     father_class + 
                     education_university_mother_dummy +
                     education_university_father_dummy + 
                     year + 
                     class5_factor + 
                     country_factor"))
  model_full<- lm(formula=formula, data=data)
  return(model_full)
}


eu_full_edu <- model_full("euftf", ESS)
immi1_full_edu <- model_full("immigration_icw_outcome", ESS_immi)
voteprog_full_edu <- model_full("vote.prog", ESS)
voterrp_full_edu <- model_full("vote.rrp", ESS)

# The full models
model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "sorting_partner_wsingle +
                     hinctnta_01 +
                     female +
                     agea10 + 
                     pdwrk + 
                     education +
                     domicil_01 +
                     ethnic_min +
                     polintr + 
                     mother_class + 
                     father_class + 
                     education_university_mother_dummy +
                     education_university_father_dummy + 
                     year + 
                     class5_factor + 
                     country_factor"))
  model_full<- lm(formula=formula, data=data)
  return(model_full)
}


eu_full <- model_full("euftf", ESS)
immi1_full <- model_full("immigration_icw_outcome", ESS_immi)
voteprog_full <- model_full("vote.prog", ESS)
voterrp_full <- model_full("vote.rrp", ESS)

## Means for the vars
eu_means <- ggemmeans(eu_full, terms=c("sorting_partner_wsingle"), rg.limit = 300000)
eu_means$outcome <- "eu"
immi_means <- ggemmeans(immi1_full, terms=c("sorting_partner_wsingle"), rg.limit = 300000)
immi_means$outcome <- "immi"
rrp_means <- ggemmeans(voterrp_full, terms=c("sorting_partner_wsingle"), rg.limit = 300000)
rrp_means$outcome <- "prog"
prog_means <- ggemmeans(voteprog_full, terms=c("sorting_partner_wsingle"), rg.limit = 300000)
prog_means$outcome <- "rrp"

eu_group_means <- ggemmeans(eu_full_edu, terms=c("education_university_dummy"), rg.limit = 300000)
immi_group_means <- ggemmeans(immi1_full_edu, terms=c("education_university_dummy"), rg.limit = 300000)
rrp_group_means <- ggemmeans(voterrp_full_edu, terms=c("education_university_dummy"), rg.limit = 300000)
prog_group_means <- ggemmeans(voteprog_full_edu, terms=c("education_university_dummy"), rg.limit = 300000)

# Add an education grouping variable
eu_means$education <- ifelse(grepl("High", eu_means$x), "Higher Educated", "Lower Educated")
immi_means$education <- ifelse(grepl("High", immi_means$x), "Higher Educated", "Lower Educated")
rrp_means$education <- ifelse(grepl("High", rrp_means$x), "Higher Educated", "Lower Educated")
prog_means$education <- ifelse(grepl("High", prog_means$x), "Higher Educated", "Lower Educated")

# Add the group means
immi_means$mean_immi <- immi_group_means$predicted[2]
immi_means$mean_immi[immi_means$education=="Lower Educated"] <- immi_group_means$predicted[1]

rrp_means$mean_rrp <- rrp_group_means$predicted[2]
rrp_means$mean_rrp[rrp_means$education=="Lower Educated"] <- rrp_group_means$predicted[1]

eu_means$mean_eu <- eu_group_means$predicted[2]
eu_means$mean_eu[eu_means$education=="Lower Educated"] <- eu_group_means$predicted[1]

prog_means$mean_prog <- prog_group_means$predicted[2]
prog_means$mean_prog[prog_means$education=="Lower Educated"] <- prog_group_means$predicted[1]


# Get the percentages of each sorting group for the labels
perc <- ESS %>% drop_na(sorting_partner_wsingle) %>% ungroup() %>% mutate(n_total = n()) %>% group_by(sorting_partner_wsingle) %>% mutate(n_group = n()) %>%
  dplyr::select(n_total, n_group) %>%  distinct(sorting_partner_wsingle, .keep_all = T)
perc$percentage <- round(perc$n_group / perc$n_total * 100, 0)


# Set the levels and labels. 
immi_means$x <- factor(immi_means$x, levels=c("Pure Low w par", "Pure Low no par", "Low ties w par", 
                                    "High ties w par", "Pure High no par", "Pure High w par"),
                  labels=c("Low, low p (37%)", "Low, no p (22%)", 
                           "Low, high p (8%)", "High, low p (9%)",
                           "High, no p (9%)", "High, high p (15%)"))
rrp_means$x <- factor(rrp_means$x, levels=c("Pure Low w par", "Pure Low no par", "Low ties w par", 
                                              "High ties w par", "Pure High no par", "Pure High w par"),
                       labels=c("Low, low p (37%)", "Low, no p (22%)", 
                                "Low, high p (8%)", "High, low p (9%)",
                                "High, no p (9%)", "High, high p (15%)"))
eu_means$x <- factor(eu_means$x, levels=c("Pure Low w par", "Pure Low no par", "Low ties w par", 
                                              "High ties w par", "Pure High no par", "Pure High w par"),
                       labels=c("Low, low p (37%)", "Low, no p (22%)", 
                                "Low, high p (8%)", "High, low p (9%)",
                                "High, no p (9%)", "High, high p (15%)"))
prog_means$x <- factor(prog_means$x, levels=c("Pure Low w par", "Pure Low no par", "Low ties w par", 
                                              "High ties w par", "Pure High no par", "Pure High w par"),
                       labels=c("Low, low p (37%)", "Low, no p (22%)", 
                                "Low, high p (8%)", "High, low p (9%)",
                                "High, no p (9%)", "High, high p (15%)"))

plot_pred_eu <- ggplot(eu_means, 
                       aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted position on EU scale",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none") +
  geom_point(aes(x=mean_eu, y=, color=education), shape=108, size=7, alpha = 0.6)

plot_pred_immi <- ggplot(immi_means, 
                         aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted position on immigration scale",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none") +
  geom_point(aes(x=mean_immi, y=x, color=education), shape=108, size=7, alpha = 0.6)

plot_pred_rrp <- ggplot(rrp_means, 
                        aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted probability to vote RRP",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none", axis.text.y = element_blank()) +
  geom_point(aes(x=mean_rrp, y=x, color=education), shape=108, size=7, alpha = 0.6)


plot_pred_prog <- ggplot(prog_means, 
                         aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted probability to vote progressive",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none", axis.text.y = element_blank()) +
  geom_point(aes(x=mean_prog, y=x, color=education), shape=108, size=7, alpha = 0.6)

plots <- grid.arrange(plot_pred_eu, plot_pred_prog, plot_pred_immi, plot_pred_rrp, nrow = 2, ncol=2)
ggsave(filename="ESS partner/plots/fig2_pred_plots_partner.png", plot=plots, dpi=600, width=9, height=5) ## Figure 2 partners








############# For parents, for each outcome, a graph with predicted probabilities
# and for the means per education group
model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "education_university_dummy +
                     hinctnta_01 +  
                      female + 
                      agea10 + 
                      pdwrk + 
                      education +
                      domicil_01 +
                      ethnic_min +
                      year + 
                      class5_factor +
                      mother_class + 
                      father_class + 
                      country_factor"))
  model_full<- lm(formula=formula, data=data)
  return(model_full)
}


eu_full_edu <- model_full("euftf", ESS)
immi1_full_edu <- model_full("immigration_icw_outcome", ESS_immi)
voteprog_full_edu <- model_full("vote.prog", ESS)
voterrp_full_edu <- model_full("vote.rrp", ESS)

# The full models
model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "sorting_parents +
                     hinctnta_01 +  
                      female + 
                      agea10 + 
                      pdwrk + 
                      education +
                      domicil_01 +
                      ethnic_min +
                      year + 
                      class5_factor +
                      mother_class + 
                      father_class + 
                      country_factor"))
  model_full<- lm(formula=formula, data=data)
  return(model_full)
}


eu_full <- model_full("euftf", ESS)
immi1_full <- model_full("immigration_icw_outcome", ESS_immi)
voteprog_full <- model_full("vote.prog", ESS)
voterrp_full <- model_full("vote.rrp", ESS)


## Means for the EU var
eu_means <- ggemmeans(eu_full, terms=c("sorting_parents"), rg.limit=29800)
eu_means$outcome <- "eu"
immi_means <- ggemmeans(immi1_full, terms=c("sorting_parents"), rg.limit=50000)
immi_means$outcome <- "immi"
rrp_means <- ggemmeans(voterrp_full, terms=c("sorting_parents"), rg.limit=50000)
rrp_means$outcome <- "prog"
prog_means <- ggemmeans(voteprog_full, terms=c("sorting_parents"), rg.limit=50000)
prog_means$outcome <- "rrp"


eu_group_means <- ggemmeans(eu_full_edu, terms=c("education_university_dummy"))
immi_group_means <- ggemmeans(immi1_full_edu, terms=c("education_university_dummy"), rg.limit=50000)
rrp_group_means <- ggemmeans(voterrp_full_edu, terms=c("education_university_dummy"), rg.limit=50000)
prog_group_means <- ggemmeans(voteprog_full_edu, terms=c("education_university_dummy"), rg.limit=50000)


# Add an education grouping variable
eu_means$education <- ifelse(grepl("High", eu_means$x), "Higher Educated", "Lower Educated")
immi_means$education <- ifelse(grepl("High", immi_means$x), "Higher Educated", "Lower Educated")
rrp_means$education <- ifelse(grepl("High", rrp_means$x), "Higher Educated", "Lower Educated")
prog_means$education <- ifelse(grepl("High", prog_means$x), "Higher Educated", "Lower Educated")


# Add the group means
immi_means$mean_immi <- immi_group_means$predicted[2]
immi_means$mean_immi[immi_means$education=="Lower Educated"] <- immi_group_means$predicted[1]

rrp_means$mean_rrp <- rrp_group_means$predicted[2]
rrp_means$mean_rrp[rrp_means$education=="Lower Educated"] <- rrp_group_means$predicted[1]

eu_means$mean_eu <- eu_group_means$predicted[2]
eu_means$mean_eu[eu_means$education=="Lower Educated"] <- eu_group_means$predicted[1]

prog_means$mean_prog <- prog_group_means$predicted[2]
prog_means$mean_prog[prog_means$education=="Lower Educated"] <- prog_group_means$predicted[1]



# Get the percentages of each sorting group for the labels
perc <- ESS %>% drop_na(sorting_parents) %>% ungroup() %>% mutate(n_total = n()) %>% group_by(sorting_parents) %>% mutate(n_group = n()) %>%
  dplyr::select(n_total, n_group) %>%  distinct(sorting_parents, .keep_all = T)
perc$percentage <- round(perc$n_group / perc$n_total * 100, 0)


# Set the levels and labels. 
immi_means$x <- factor(immi_means$x, levels=c("Pure Low", "Low weak ties", "Low strong ties", 
                                              "High strong ties", "High weak ties", "Pure High"),
                       labels=c("Low, 0 (60%)", "Low, 1 (5%)", "Low, 2 (1%)", "High, 2 (21%)", "High, 1 (8%)", "High, 0 (5%)"))
rrp_means$x <- factor(rrp_means$x, levels=c("Pure Low", "Low weak ties", "Low strong ties", 
                                              "High strong ties", "High weak ties", "Pure High"),
                       labels=c("Low, 0 (60%)", "Low, 1 (5%)", "Low, 2 (1%)", "High, 2 (21%)", "High, 1 (8%)", "High, 0 (5%)"))
eu_means$x <- factor(eu_means$x, levels=c("Pure Low", "Low weak ties", "Low strong ties", 
                                              "High strong ties", "High weak ties", "Pure High"),
                       labels=c("Low, 0 (60%)", "Low, 1 (5%)", "Low, 2 (1%)", "High, 2 (21%)", "High 1, (8%)", "High, 0 (5%)"))
prog_means$x <- factor(prog_means$x, levels=c("Pure Low", "Low weak ties", "Low strong ties", 
                                              "High strong ties", "High weak ties", "Pure High"),
                       labels=c("Low, 0 (60%)", "Low, 1 (5%)", "Low, 2 (1%)", "High, 2 (21%)", "High, 1 (8%)", "High, 0 (5%)"))


plot_pred_eu <- ggplot(eu_means, 
                       aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted position on EU scale",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none") +
  geom_point(aes(x=mean_eu, y=, color=education), shape=108, size=7, alpha = 0.6)

plot_pred_immi <- ggplot(immi_means, 
                         aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted position on immigration scale",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none") +
  geom_point(aes(x=mean_immi, y=x, color=education), shape=108, size=7, alpha = 0.6)

plot_pred_rrp <- ggplot(rrp_means, 
                        aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted probability to vote RRP",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none", axis.text.y = element_blank()) +
  geom_point(aes(x=mean_rrp, y=x, color=education), shape=108, size=7, alpha = 0.6)


plot_pred_prog <- ggplot(prog_means, 
                         aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted probability to vote progressive",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) + theme(legend.position = "none", axis.text.y = element_blank()) +
  geom_point(aes(x=mean_prog, y=x, color=education), shape=108, size=7, alpha = 0.6)



plots <- grid.arrange(plot_pred_eu, plot_pred_prog, plot_pred_immi, plot_pred_rrp, nrow = 2, ncol=2)
ggsave(filename="ESS partner/plots/fig2_pred_plots_parents.png", plot=plots, dpi=600, width=9, height=5) ### Figure 2 parents













############# Robustness analysis: results by country ##########################
#### Partner country
# A model where we interact the treatment variable with parents
country_interact <- lm(immigration_icw_outcome~ sorting_partner_wsingle * country_factor +
                         hinctnta_01 +
                         female +
                         agea10 + 
                         pdwrk + 
                         education +
                         domicil_01 +
                         ethnic_min +
                         polintr + 
                         mother_class + 
                         father_class + 
                         education_university_mother_dummy +
                         education_university_father_dummy + 
                         year + 
                         class5_factor 
                       , data=ESS_immi)



# And a model to get the country means on immigration with our other variables included
country_interact_edu <- lm(immigration_icw_outcome~ education_university_dummy * country_factor +
                             hinctnta_01 +
                             female +
                             agea10 + 
                             pdwrk + 
                             education +
                             domicil_01 +
                             polintr + 
                             ethnic_min +
                             mother_class + 
                             father_class + 
                             education_university_mother_dummy +
                             education_university_father_dummy + 
                             year + 
                             class5_factor
                           , data=ESS_immi)




immi_means <- ggemmeans(country_interact, terms=c("sorting_partner_wsingle", "country_factor"), rg.limit=50000) #by country, predicted value for the treatment
immi_means$outcome <- "immi"
immi_means$education <- ifelse(grepl("High", immi_means$x), "Higher Educated", "Lower Educated") # add education grouping variable

# calculate the group means
immi_group_means <- ggemmeans(country_interact_edu, terms=c("education_university_dummy", "country_factor"), rg.limit=50000)
immi_group_means <- as_tibble(immi_group_means)
immi_group_means$education <- ifelse(immi_group_means$x==1, "Higher Educated", "Lower Educated")
immi_group_means <- immi_group_means %>% rename(mean_immi=predicted) %>% dplyr::select(mean_immi, group, education)

# add the mean level of education to the country data
immi_means <- full_join(immi_means, immi_group_means)

immi_means$x <- factor(immi_means$x, levels=c("Pure Low w par", "Pure Low no par", "Low ties w par", 
                                              "High ties w par", "Pure High no par", "Pure High w par"),
                       labels=c("Low, low p", "Low, no p", 
                                "Low, high p", "High, low p",
                                "High, no p", "High, high p"))


partner_countries <- ggplot(immi_means, aes(x=predicted, y=x, xmin=conf.low, xmax=conf.high, color=education, shape=education)) + 
  geom_pointrange(size=0.8) + 
  labs(y="", x="Predicted position on immigration scale",
       color="Education or respondents:", shape="Education or respondents:",
       title="") +
  scale_color_grey(end=0.6) +
  geom_point(aes(x=mean_immi, y=x, color=education), shape=108, size=7, alpha = 0.6) +
  facet_wrap(~group, scales = "free") +
  theme(legend.position = "none", 
        strip.text = element_text(colour = 'black'),
        strip.background =element_rect(fill="lightgrey")
  )

### Plot A1
ggsave(plot=partner_countries, "ESS partner/plots/figa1_partner_countries.png", dpi = 600, width=10, height=7, device="png")




























############# Robustness analysis: other immigration variables ##########################


##### Main table with all other immigration outcomes in the ESS for partners.
model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "sorting_partner_wsingle +
                     hinctnta_01 +
                     female +
                     agea10 + 
                     pdwrk + 
                     education +
                     domicil_01 +
                     ethnic_min +
                     year + 
                     class5_factor + 
                     country_factor"))
  model_full<- lm(formula=formula, data=data)
  return(model_full)
}


## Simple Additive index of six different immigration items: 
# imsmetn   imdfetn   impcntr   imbgeco   imueclt   imwbcnt

## Reverse the polarity of the items where higher values are not more progressive
ESS$imsmetn_rev <- max(ESS$imsmetn,na.rm=T) - ESS$imsmetn
ESS$imdfetn_rev <- max(ESS$imdfetn,na.rm=T) - ESS$imdfetn
ESS$impcntr_rev <- max(ESS$impcntr,na.rm=T) - ESS$impcntr



# Country's cultural life undermined or enriched by immigrants
immi1_full <- model_full("imueclt", ESS)

#Immigration bad or good for country's economy
immi2_full <- model_full("imbgeco", ESS)

immi3_full <- model_full("imwbcnt", ESS)

immi4_full <- model_full("impcntr_rev", ESS)

immi5_full <- model_full("imsmetn_rev", ESS)

#Allow many/few immigrants of different race/ethnic group from majority
immi6_full <- model_full("imdfetn_rev", ESS)

#cavaille
immi7_full <- model_full("immi_pro_cavaille", ESS)



### And running the models for parents
model_full <- function(outcomevar, data){
  formula <- formula(paste0(outcomevar, "~", "sorting_parents + 
                      hinctnta_01 +  
                      female + 
                      agea10 + 
                      pdwrk + 
                      education +
                      domicil_01 +
                      ethnic_min +
                      year + 
                      class5_factor +
                      mother_class + 
                      father_class + 
                      country_factor"))
  model_full <- lm(formula=formula, data=data)
  return(model_full)
}



# Country's cultural life undermined or enriched by immigrants
immi1_full_parents <- model_full("imueclt", ESS)

#Immigration bad or good for country's economy
immi2_full_parents <- model_full("imbgeco", ESS)

immi3_full_parents <- model_full("imwbcnt", ESS)

immi4_full_parents <- model_full("impcntr_rev", ESS)

immi5_full_parents <- model_full("imsmetn_rev", ESS)

#Allow many/few immigrants of different race/ethnic group from majority
immi6_full_parents <- model_full("imdfetn_rev", ESS)

#cavaille
immi7_full_parents <- model_full("immi_pro_cavaille", ESS)





##### Main table with all immigration outcomes for parents
model <- texreg(l=list(immi1_full, immi2_full, immi3_full, immi4_full, immi5_full, immi6_full, immi7_full, 
                       immi1_full_parents, immi2_full_parents, immi3_full_parents, immi4_full_parents, immi5_full_parents, immi6_full_parents, immi7_full_parents),
                stars=c(0.05, 0.01, 0.001),use.packages=F,
                digits=3, include.ci = FALSE, booktabs=T,
                custom.model.names = c("imueclt", "imbgeco", "imwbcnt", "impcntr", "imsmetn", "imdfetn", "CM", 
                                       "imueclt", "imbgeco", "imwbcnt", "impcntr", "imsmetn", "imdfetn", "CM"),
                omit.coef=c("(year)|(employed)|(urban)|(migration)|(agea10)|(father)|(mother)|(country)|(hinctnta_01)|(total)|(female)|(pdwrk)|(ducation)|(domicil_01)|(class5)|(mother)|(father)|(country)|(ethnic_min)"), 
                custom.coef.names=c("Intercept", "High res no par", "High res Low par",
                                    "Low res High par", "Low res no par", 
                                    "Low res Low par",
                                    "H edu One L parent", "High edu L parents", "Low edu H parents", 
                                    "L edu One L parent", "Low edu L parents"),
                reorder.coef=c(2,3,4,5,6,7,8,9,10,11,1),
                caption.above=T,
                fontsize="scriptsize",
                caption="ESS: Partner and parental network effects on other immigration outcomes",
                threeparttable = TRUE,
                label="tab:ess_immigration_robust",
                float.pos="htpb!",
                custom.gof.rows=list("Controls"= c("Yes", "Yes", "Yes","Yes", "Yes", "Yes", "Yes","Yes", "Yes", "Yes","Yes", "Yes","Yes","Yes"),
                                     "Year FE"= c("Yes", "Yes", "Yes","Yes","Yes", "Yes", "Yes","Yes", "Yes", "Yes","Yes", "Yes","Yes","Yes"),
                                     "Country FE"= c("Yes", "Yes", "Yes","Yes","Yes", "Yes", "Yes","Yes", "Yes", "Yes","Yes", "Yes","Yes","Yes"),
                                     "Class FE"= c("Yes", "Yes", "Yes","Yes","Yes", "Yes", "Yes","Yes", "Yes", "Yes","Yes", "Yes","Yes","Yes")),
                custom.note="\\item %stars.")

## Table A7
print(model, file = "ESS partner/tables/taba7_ess_immigration_robust.tex")

         







